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Abstract 

This  paper  is  concerned  with  the  simultaneous  flow  of  liquid  water  and  gaseous  reactants  in  mini-channels  of  a  proton  exchange  membrane 
(PEM)  fuel  cell.  Envisaging  the  mini-channels  as  structured  and  ordered  porous  media,  we  develop  a  continuum  model  of  two-phase  channel 
flow  based  on  two-phase  Darcy’s  law  and  the  M2  formalism,  which  allow  estimate  of  the  parameters  key  to  fuel  cell  operation  such  as  overall 
pressure  drop  and  liquid  saturation  profiles  along  the  axial  flow  direction.  Analytical  solutions  of  liquid  water  saturation  and  species  concentrations 
along  the  channel  are  derived  to  explore  the  dependences  of  these  physical  variables  vital  to  cell  performance  on  operating  parameters  such  as 
flow  stoichiometric  ratio  and  relative  humility.  The  two-phase  channel  model  is  further  implemented  for  three-dimensional  numerical  simulations 
of  two-phase,  multi-component  transport  in  a  single  fuel-cell  channel.  Three  issues  critical  to  optimizing  channel  design  and  mitigating  channel 
flooding  in  PEM  fuel  cells  are  fully  discussed:  liquid  water  buildup  towards  the  fuel  cell  outlet,  saturation  spike  in  the  vicinity  of  flow  cross-sectional 
heterogeneity,  and  two-phase  pressure  drop.  Both  the  two-phase  model  and  analytical  solutions  presented  in  this  paper  may  be  applicable  to  more 
general  two-phase  flow  phenomena  through  mini-  and  micro-channels. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Fuel  cells,  converting  chemical  energy  of  fuels  directly 
into  electricity,  have  become  an  integral  part  of  alternative 
energy  and  energy  efficiency.  Their  noteworthy  features,  high- 
energy  conversion  efficiency  and  zero  emission,  meet  the 
critical  demands  of  a  rapidly  growing  society  [1,2].  Among 
all  types  of  fuel  cells,  the  proton  exchange  membrane  (PEM) 
fuel  cell,  also  called  polymer  electrolyte  fuel  cell  (PEFC), 
has  reached  center  stage,  particularly  for  mobile  and  portable 
applications  [3,4].  Besides  their  high-power  capability,  PEM 
fuel  cells  work  at  low  temperatures,  produce  only  water  as 
byproduct,  and  can  be  compactly  assembled,  making  them 
one  of  the  leading  candidates  for  the  next  generation  power 
generator. 

A  typical  PEFC  consists  of  bipolar  plates,  gas  channels,  gas 
diffusion  layers  (GDLs),  and  a  proton-conductive  membrane 
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with  platinum  catalyst  coated  on  each  side,  called  the  membrane 
electrode  assembly  (MEA),  as  shown  in  Fig.  1.  Gas  channels 
are  grooved  in  graphite  or  metal  plates,  where  injected  reactant 
streams  are  distributed  for  electrochemical  reactions.  The  GDLs, 
usually  coated  with  micro-porous  layers  (MPLs),  play  an  impor¬ 
tant  role  in  electronic  connection  between  the  bipolar  plate  and 
the  electrode  and  provide  a  passage  for  reactant  transport  and 
heat/water  removal.  Protons  are  produced  from  hydrogen  oxida¬ 
tion  in  the  anode  catalyst  layer,  and  pass  through  the  membrane, 
carrying  water  molecules  via  electro-osmotic  drag,  to  the  cath¬ 
ode  catalyst  layer  where  the  oxygen  reduction  reaction  (ORR) 
occurs  with  water  as  byproduct. 

Water  management  is  a  central  issue  in  PEFC  technology 
because  while  water  is  essential  for  membrane  ionic  conduc¬ 
tivity,  excess  liquid  water  leads  to  flooding  of  catalyst  layers 
and  GDLs  [5-7]  as  well  as  channel  clogging  [8,9].  Given  low- 
operating  temperatures  during  normal  startup  (25  °C)  and  hence 
low-saturation  pressures,  two-phase  phenomena  are  unavoid¬ 
able  in  automotive  fuel  cells.  Two-phase  transport  in  a  fuel  cell 
consists  of  three  sub-problems:  catalyst  layer  flooding,  GDL 
flooding,  and  two-phase  flow  in  channels.  To  date,  most  of  the 
efforts  in  two-phase  modeling  were  devoted  to  the  former  two 
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Fig.  1.  Schematic  of  a  PEM  fuel  cell. 


[5,10-18],  and  few  models  for  channel  flooding  have  existed  in 
the  literature. 

In  practice,  however,  two-phase  flow  and  transport  in  chan¬ 
nels  is  of  paramount  importance  for  fuel  cell  operation, 
particularly  in  the  most  energy-efficient  regime  involving  high¬ 
cell  voltage  and  therefore  low-current  density.  In  this  regime, 
flow  rates  of  hydrogen  and  air  through  channels  are  exceed¬ 
ingly  low,  resulting  in  substantial  liquid  water  accumulation. 
Fig.  2  shows  experimental  data  of  cell  voltage  and  cathode-side 
pressure  drop  as  functions  of  air  stoichiometry.  It  can  be  seen 
that  channel  flooding  at  low  stoichiometry  (§  =  2)  results  in  cell 
voltage  drop  of  as  much  as  120  mV,  which  would  negate  a  volt¬ 
age  gain  of  45  mV  targeted  by  fourfold  improvement  in  catalyst 
activity.  Moreover,  the  drastic  fluctuation  in  cell  voltage  at  low 
stoichiometry  sets  up  voltage  cycling  at  high  potentials,  which 
has  severe  impact  on  PEM  fuel  cell  durability.  Other  major  con¬ 
cerns  of  channel  flooding  include:  (1)  fuel  starvation  which  leads 
to  carbon  corrosion  in  the  cathode  catalyst  layer;  (2)  oxygen  star¬ 
vation  which  leads  to  hydrogen  evolution  on  the  cathode  and 
furthermore  H2/O2  combustion  and  hot  spot  formation;  (3)  flow 
maldistribution  among  parallel  channels,  leading  to  operational 


instability  and  efficiency  losses;  (4)  increase  of  the  mass  trans¬ 
port  loss  at  high-current  densities.  Thus,  an  urgent  need  in  PEFC 
modeling  has  been  to  develop  a  two-phase  flow  and  flooding 
model  for  fuel  cell  channels,  allowing  the  liquid  water  satura¬ 
tion  to  be  predicted  at  levels  as  high  as  15-20%.  Developing 
such  a  preliminary  model  is  the  objective  of  this  work. 

In  this  work,  we  envision  a  structural  and  flow  analogy 
between  PEM  fuel  cell  channels  and  random  porous  media. 
Following  the  analogy,  we  then  apply  two-phase  flow  theory 
in  porous  media  to  describe  the  two-phase,  multi-component 
transport  in  channels  of  a  fuel  cell.  We  show  that  such  a 
phenomenological  model  is  able  to  capture  pressure,  liquid  sat¬ 
uration  and  reactant  composition  distributions  along  the  flow, 
thereby  sufficing  to  address  such  questions  as  channel  flooding, 
water  trapping  in  channels  of  heterogeneous  geometry,  two- 
phase  flow  maldistribution  in  multiple,  parallel  channels,  and 
the  flowfield  effect  on  liquid  water  drainage.  In  addition,  we 
present  analytical  solutions  as  to  the  saturation  distribution  along 
the  channel  and  analysis  on  the  heterogeneity  impact  and  two- 
phase  pressure  factor  for  PEFC  channels  for  the  first  time  as 
well  as  three-dimensional  numerical  simulations  to  explore  the 
effects  of  key  parameters  such  as  relative  humidity  (RH)  and 
air  stoichiometric  ratio  on  two-phase  flow  in  channels.  Finally, 
experimental  validation  of  the  predicted  two-phase  pressure 
drop  is  reported. 

2.  Analogy  between  mini-channels  and  random  porous 
media 

PEM  fuel  cell  channels  are  characterized  by  parallel  or 
serpentine  channels  with  square,  rectangular  or  trapezoid  cross- 
sections  and  channel  dimensions  ranging  from  0.1  to  1mm 
(see  Figs.  3  and  4)  [2,19].  Different  configurations  have  been 
described  and  explored  recently,  such  as  serpentine,  parallel, 
and  interdigitated  flow  fields  [19-23].  Fig.  4  shows  a  typical 
serpentine  parallel  flow  field  in  an  industrial  PEFC  [24].  Visual¬ 
ization  studies  have  revealed  and  studied  two-phase  phenomena 
in  fuel  cell  channels  experimentally  [8,9,25],  as  illustrated  in 
Fig.  5.  Channels  of  similar  structure  are  also  encountered  in 
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Fig.  2.  Cell  voltage  and  cathode-side  pressure  drop  measured  under  various  air  stoichiometric  ratios  (symbolized  by  £  in  the  figure)  in  a  14  cm2  fuel  cell  for  0.2  A  cm-2 , 
150  kPa  and  80  °C  with  the  dew  point  of  both  H2  and  air  at  70  °C  [54] .  As  the  air  stoichiometry  decreases,  the  pressure  drop  does  not  decrease  proportionally,  indicating 
the  onset  of  severe  water  flooding  in  cathode  channels.  Simultaneously,  the  cell  voltage  drops  by  as  much  as  120  mV  and  fluctuates  more. 
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(C)  <d) 

Fig.  3.  Typical  flow  fields  of  PEFCs.  (a)  Parallel  flow  field;  (b)  serpentine  flow  field;  (c).  pin  flow  field;  (d)  interdigitated  flow  field. 


Outlet  manifold 


Fig.  4.  A  flow  field  of  an  industrial  PEFC  [24]. 


miniature  heat  pipes  and  micro-heat  exchangers.  In  these  appli¬ 
cations,  evaporation  and  condensation  take  place  to  utilize  the 
efficient  latent  heat  transfer,  therefore  two-phase  flow  is  ubiqui¬ 
tous.  Previous  studies  in  these  fields  include  those  of  Stubos  et 
al.  [26],  Wang  et  al.  [27],  Tio  et  al.  [28],  Imke  [29],  Henning  et 
al.  [30],  and  Sugumar  and  Tio  [31]. 

Another  field  worthy  mentioning  is  two-phase  flow  in  geo¬ 
logical  media  or  petroleum  reservoirs.  Flow  beds  in  these 
applications  feature  random  solid  particles  with  two-phase  flow 
in  pores.  Some  porous  formations  feature  a  pore  size  range  of 
~lmm,  thus  strongly  resembling  the  fuel  cell  channels.  The 
structural  similarity  between  the  two,  with  the  former  being 
random  porous  media  and  the  latter  regular,  ordered  pores,  was 
recognized  by  Chaouche  et  al.  [32],  Or  and  Tuller  [33],  and  Li  et 
al.  [34] .  In  addition,  there  exists  a  flow  analogy  in  that  flow  in  fuel 
cell  channels  is  well  within  the  laminar  regime  and  exhibits  a  lin¬ 
ear  relationship  between  the  pressure  drop  and  velocity,  the  same 
relation  as  expressed  by  Darcy’s  law  for  flow  in  porous  media. 
Indeed,  the  proportionality  constant  called  hydraulic  conduc¬ 
tance  in  laminar  flow  through  channels  is  physically  equivalent 
to  permeability  of  a  porous  medium.  Such  a  flow  analogy  was 
recognized  by  Chen  [35],  Wong  et  al.  [36],  Wang  et  al.  [27],  and 


Fig.  5.  Visualization  of  two-phase  flow  in  PEM  fuel  cell  channels  [55]. 
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Tio  et  al.  [28],  among  others.  In  addition,  the  geometrical  anal¬ 
ogy  between  a  porous  medium  and  a  capillary  is  traditionally 
employed  to  understand  and  develop  flow  theories  in  irregular 
porous  media,  in  which  disordered  porous  media  with  tortu¬ 
ous  pore  channels  are  simplified  as  a  bundle  of  straight  capillary 
tubes.  Envisaging  the  geometrical  and  flow  analogy  between  fuel 
cell  channels  and  random  porous  media,  we  seek  to  develop  a 
macroscopic  model  for  two-phase  flow  in  fuel  cell  channels  by 
using  the  two-phase  Darcy’s  flow  theory,  as  detailed  in  the  next 
section.  Similar  approaches  were  successfully  implemented  by 
Wang  et  al.  [27],  Tio  et  al.  [28],  and  Imke  [29]  for  miniature  heat 
pipes  and  micro-heat  exchangers. 

3.  Physical  model 

Fig.  1  schematically  shows  the  geometry  of  a  PEFC  and  its 
components.  The  focus  is  placed  on  a  single  straight  channel  in 
order  to  develop  a  channel  two-phase  flow  model  resolving  axial 
variations  in  pressure,  liquid  saturation  and  other  quantities.  This 
is  schematically  shown  in  Fig.  6(a)  for  the  cathode  side.  Partially 
or  fully  humidified  air  is  fed  in  the  inlet,  and  liquid  water  pro¬ 
duced  from  ORR  is  injected  into  the  channel  from  the  sidewall 
facing  the  GDF.  The  channels  can  generally  be  treated  as  porous 
media  with  or  without  porous  inserts.  Without  porous  inserts,  the 
pore  size  in  the  porous  medium  is  exactly  equal  to  the  channel 
cross-section  dimension  with  the  pore  tortuosity  and  porosity 
equal  to  unity.  This  general  porous  approach  is  illustrated  in 
Fig.  6(b)  and  the  Darcy’s  law  is  applied  to  describe  the  two-phase 
flow  through  the  channels,  provided  that  the  Reynolds  number 
based  on  the  pore  size  is  much  smaller  than  2000  to  ensure  lami¬ 
nar  flow.  This  condition  is  commonly  met  in  fuel  cell  operation. 
Also,  the  capillary  number  denoting  the  ratio  of  viscous  force 
to  surface  tension  is  unimportant  for  channel  two-phase  flow  in 
the  presence  of  large  axial  velocity  in  gas  channels. 

3.1.  Conservation  of  mass  and  momentum 


for  two-phase  flow  through  porous  media  [37],  the  mass  and 
momentum  conservations  of  the  two-phase  mixture  can  be  writ¬ 
ten  as  follows: 

dp 

Continuity  equation  :  s - f  V  •  ( pu )  =  0  (1) 

3 1 

g 

Momentum  conservation  :  pu  = - (VP  —  yppg)  (2) 

v 

where  the  physical  properties  of  the  two-phase  mixture  are 
defined  as 


p  =  spl  +  (l-  s)pg 

pU  U\P\  “l-  UgPg 


V  = 


^rl  &rg\ 

VI  Vg) 
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Here,  the  liquid  water  saturation,  s,  is  defined  as  the  volume 
fraction  of  open  pores  occupied  by  liquid  water,  which  is  a  key 
parameter  characterizing  the  two-phase  flow.  Fiquid  saturation 
is  identical  to  liquid  volume  fraction  in  channel  two-phase  flow. 
The  liquid  saturation  can  be  determined  from  the  following  rela¬ 
tion  with  the  mixture  water  concentration,  CH2°,  after  the  latter 
is  solved  from  the  water  transport  equation: 

CHl°  -  Csa, 

S  Pl/MW  -  Csat  } 


The  density  correction  factor  in  Eq.  (2),  yp,  is  a  sole  function 
of  phase  saturations  and  thus  can  be  regarded  as  a  property  of 
the  two-phase  mixture: 


_  Pih  +  Pg^g 
Yp  sp\  +  (1  -  s)pg 

where  the  phase  mobility,  A&,  is  defined  as 


(7) 
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and  Ag  =  1  —  X\ 


(8) 


In  the  channel  pores,  we  consider  the  gas  and  liquid  phases  as 
constituents  of  a  two-phase  mixture  with  distinct  phase  veloci¬ 
ties.  Following  the  well-known  multiphase  mixture  formulation 
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Fig.  6.  Schematic  of  two-phase  flow  through  a  single  channel:  (a)  an  open  chan¬ 
nel;  (b)  a  channel  with  porous  inserts.  The  black  shadowed  regions  on  the  upper 
or  lower  wall  represent  the  liquid  phase. 


The  relationship  between  the  phase  mobility,  A&,  and  liquid 
water  saturation,  s,  is  plotted  in  Fig.  7.  It  can  be  seen  that  the 


Fig.  7.  Phase  mobility  vs.  liquid  water  saturation  (; n k  =  4). 
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liquid  mobility  increases  with  liquid  saturation.  The  relative  per¬ 
meability,  krg,  is  defined  as  the  ratio  of  the  intrinsic  permeability 
of  liquid/gas  phase  to  the  total  permeability  of  a  porous  medium. 
Physically,  it  describes  the  extent  to  which  one  fluid  is  hindered 
by  others  in  pore  spaces  and  hence  can  be  formulated  as  a  func¬ 
tion  of  liquid  saturation.  A  common  set  of  functions  have  been 
used  in  the  literature: 


kri  = 


{  S  ~  % 


nk 


and 


(9) 


Beming  and  Djilali  [38]  set  ng  equal  to  unity  in  their  work, 
a  linear  relationship  also  called  X  curve,  while  some  used  the 
value  of  3  or  4  [6, 1 3 , 1 6, 1 7] .  Note  that  the  larger  ng  is,  the  greater 
the  degree  to  which  the  liquid  water  affects  the  gas  flow.  In  the 
present  work,  we  use  ng  =  4  following  Corey  [39]. 

The  irreducible  residual  saturation,  S{r,  is  physically  the  liquid 
fraction  bound  inside  the  capillary  channel  by  surface  tension 
force  and  hence  cannot  be  removed  by  drainage  but  only  by  evap¬ 
oration.  The  value  of  s jr  is  determined  by  the  balance  between 
surface  tension  force,  a,  and  gravitational  force,  and  is  suggested 
by  Saez  and  Carbonell  [40]  for  mini-channels  in  the  follow 
function  of  Eotvos  number,  Eo: 


1 

20  +  0.9E  o 


(10) 


The  Eotvos  number  is  defined  by: 


Ed  = 


Pigdj 

a 


(ID 


where  <4  is  the  hydraulic  diameter  of  the  channel. 

Once  kvg  is  available,  the  following  Darcy’s  law  is  applied  to 
calculate  the  velocity  for  each  phase: 

PkMk  =  -krk—(yPk  -  Pkg)  (12) 

Vk 

In  addition,  a  dimensionless  number,  Bo,  is  defined  as  the 
ratio  of  the  gravitational  force  to  the  surface  tension,  Bo  = 
(g  Ap/cr)d^  where  Ap,  <4  and  a  are  the  density  difference 
between  phases,  the  channel  hydraulic  diameter  and  surface 
tension,  respectively.  Bo  is  typically  small  in  fuel  cells  (<0.1); 
therefore,  the  gravitational  term  in  Eq.  (12)  can  be  neglected. 

Cathode  gas  channels  commonly  feature  a  constant  flow 
as  Wang  and  Wang  demonstrated  that  the  flow  variations 
due  to  ORR  and  phase  change  is  less  than  4%  [45].  Under 
this  assumption,  the  porous  medium  permeability,  X,  can  be 
determined  from  the  hydraulic  conductance  of  laminar  flow 
through  a  flow  channel  of  various  cross-sections.  It  is  usually 
obtained,  under  single-phase  condition,  through  the  well-known 
Hagen-Poiseuille  equation: 


K  = 


(13) 


Table  1 

Flow  shape  factors  for  typical  channel  cross-sections  in  PEM  fuel  cells  [27] 
Shape  c 


and  c  is  a  factor  determined  by  the  cross-section  shape,  which 
is  also  called  flow  shape  factor.  For  convenience,  some  c  values 
have  been  listed  in  Table  1  for  several  typical  cross-sections  of 
fuel  cell  channels. 


3.2.  Conservation  of  species 


The  primary  species  in  PEFC  channels  include  hydrogen, 
water,  and  nitrogen  in  the  anode,  and  oxygen,  nitrogen,  and  water 
in  the  cathode  side.  A  general  form  of  species  transport  equation 
for  both  single-  and  two-phase  mixtures  can  be  expressed  by 
[37]: 


Species  conservation  : 


+  V  •  (ycuCk) 


=  V  •  (Dg  effVCg)  -  V  • 


(15) 


where  yc  is  called  the  convection  correction  factor  to  correct  the 
convective  transport  of  the  two-phase  mixture  due  to  difference 
between  phase  velocities.  Assuming  that  other  species,  such  as 
oxygen,  hydrogen,  and  nitrogen,  have  a  negligible  solubility  in 
the  liquid  water,  yc  can  be  expressed  as 


Yc  = 


P  f  h 
CH 2°  yM1^0 

Mg 

pg(i  -sf 


for  water 

for  other  species 

(16) 


where  the  hydraulic  diameter  of  the  channel  is  given  by 


4  =  4x 


cross-section  area 


(14) 


In  Eq.  (15),  the  first  term  on  the  right  side  describes  the  effect 
of  the  molecular  diffusive  transport.  As  more  than  two  species 
are  present  in  a  fuel  cell,  the  multi-component  Stefan-Maxwell 
diffusion  is  usually  applied.  To  simplify  the  model  expression 


channel  perimeter 
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and  later  numerical  implementation,  we  use  the  Fick’s  law  to 
describe  the  molecular  diffusion  in  the  gaseous  phase.  In  addi¬ 
tion,  Eq.  (15)  also  neglects  diffusion  of  hydrogen,  oxygen,  and 
nitrogen  in  liquid  water  due  to  exceedingly  low  solubility.  In  the 
gas  phase,  the  species  diffusion  coefficient  in  Eq.  (15)  is  given 
by  [41]: 


rr„ 


(17) 


In  a  two-phase  region,  the  effective  gaseous  diffusion  coeffi¬ 
cient  is  modified  via  the  Bruggeman  correlation: 

Dk’eff  =  (1  -s)TiDkg  (18) 


In  addition  to  the  convective  and  diffusive  transport,  capillary 
action  may  take  place  in  porous  media,  which  drives  liquid  water 
from  a  higher  liquid  saturation  region  to  a  lower  one.  The  water 
flux  driven  by  the  capillary  force  is  called  the  capillary  diffusion 
flux,  ji,  which  is  included  in  the  second  term  of  Eq.  (15).  The 
capillary  force  is  mathematically  described  by  capillary  pres¬ 
sure,  pc,  which  is  defined  as  the  pressure  difference  between 
gas  and  liquid  phases,  i.e.  Pg—p\.  The  capillary  pressure  is  a 
function  of  surface  properties  such  as  surface  tension,  a,  and 
contact  angle,  0C,  porous  medium  properties  such  as  porosity,  s, 
and  permeability,  K,  and  the  liquid  saturation,  s.  The  capillary 
diffusion  flux,  j\,  and  capillary  pressure,  pc ,  can  be  calculated 
by 

jl  =  —K[VPc  +  (Pl~P g)|]  (19) 

V 


where 
Pc  =  cr 


(  E  \  1/2 

cos(6c)[—)  J(S ) 


Here,  J(s)  is  the  Leverett  function,  an  empirical  relation  that  is 
generally  adopted  for  both  hydrophobic  and  hydrophilic  porous 
media  [7],  and  we  extend  it  to  structured  porous  medium  repre¬ 
sentation  of  gas  channels: 


A  general  thermal  model  for  two-phase  flow  can  be  written 
as 


Energy  conservation  : 
=  V  •  (keffVT)  +  Sfg 


dpcvT 

dt 


+  V  •  ( yTpcvuT ) 


(22) 


Once  again,  the  advection  term  is  modified  by  a  correction 
factor,  yt,  and  is  given  by 

yT  —  +  7,gCP;g)  ^3) 

sp\cv,\  +  (1  -  s)pgCp,g 

The  source  term  arises  from  latent  heat  release  or  absorption 
due  to  water  condensation/evaporation  and  can  be  written  as  [  1 6] 

Sfg  =  hfgriifg  (24) 

where  hfg  and  mfg  are  the  latent  heat  of  vapor-liquid  phase 
change  and  the  phase  change  rate,  respectively.  The  latter  is 
readily  calculated  from  the  liquid  continuity  equation,  namely 

Wfg  =  p\  —  +  v  •  (A«i)  (25) 

where  the  liquid-phase  velocity  in  the  M 2  model  is  computed 
from  either  Eq.  (12)  or  the  following  equation: 

pm  =  j\  +  h  pu  (26) 


3.4.  Boundary  conditions 

Eqs.  (1),  (2),  (15)  and  (22)  form  a  complete  set  of  govern¬ 
ing  equations  with  unknowns  of  u,  P,  T,  and  Ck  to  describe 
non-isothermal,  two-phase  flow  in  fuel  cell  channels.  Proper 
boundary  conditions  are  needed  to  close  the  mathematical  sys¬ 
tem.  The  temperature  distribution  is  usually  determined  by  the 


J(s)  = 


1.417(1  -s)  -  2.120(1  -  s )2  +  1.263(1  -  s)3,  for  <9C  <  90° 
1.417s  —  2.120s2  +  1.263s3,  for  9C  >  90° 


(20) 


Usually,  the  GDL  materials  are  made  hydrophobic  through 
adding  PTFE  to  facilitate  the  water  removal  from  the 
channel-GDL  interface  [42].  Other  channel  walls  are  usually 
hydrophilic,  allowing  wicking  of  liquid  water  away  from  the 
GDL  surface. 

3.3.  Conservation  of  energy 

Temperature  distribution  plays  a  vital  role  in  channel  two- 
phase  flow  due  to  the  exponential  relationship  between  the 
saturation  pressure,  pSSL t,  and  temperature.  Springer  et  al.  [43] 
correlated  ps at  with  temperature: 


cooling  channel  design  in  PEFCs,  as  detailed  previously  by 
Wang  and  Wang  [24].  Other  boundary  conditions  are  specified 
as  follows. 


3.4.1.  Flow  inlet  boundaries 

The  cathode  inlet  flow  is  commonly  set  partially  or  fully  satu¬ 
rated  without  liquid  water.  Therefore,  the  inlet  mixture  velocity 
is  the  one  in  the  gas  phase,  um^,  determined  by  the  stoichio¬ 
metric  flow  ratio,  §c,  and  a  reference/average  current  density, 


,  _  C®2u\n  gAxz 

sc  — 


7av^mem/  4  F 


(27) 


log10  Psat  =  -2.1794  +  0.02953(7  -  273.15) 

-9.1837  x  10-5(T  —  273. 15)2 
+  1.4454  x  l(T7(r  -  273. 15)3  (21) 


where  Axz  and  Amem  are  the  areas  of  flow  cross-sectional  and 
the  membrane,  respectively.  The  inlet  molar  concentrations  are 
determined  by  the  inlet  pressure  and  humidity  according  to  the 
ideal  gas  law. 
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3.4.2.  Outlet  boundaries 

Fully  developed  conditions  are  applied: 

du  ,  dCk 

—  =0  and  - =  0 

3  n  3  n 


(28) 


3.4.3.  Solid  wall 

Slip  and  impermeable  velocity  condition  and  no- flux  condi¬ 
tion  are  applied: 


u  •  n  =  0  and 


3  Ck 
3  n 


0 


(29) 


3.4.4.  GDL  surface 

As  only  the  flow  channel  is  considered,  during  operation 
product  water  is  injected  into  the  channel  at  the  GDL-channel 
interface  and  simultaneously  a  proportional  amount  of  oxygen 
is  transported  out  of  channel  for  the  ORR  in  the  catalyst  layer. 
Given  a  distribution  of  current  density,  I(x,y),  one  can  calcu¬ 
late  the  water  production  and  oxygen  consumption  rates  via  the 
Faraday’s  law: 


h2o  _  (1  +  2a) I 
IF 


and 


_  _J_ 

4  F 


(30) 


where  a  is  the  net  water  transport  coefficient  which  describes 
the  total  water  crossover  from  the  anode  to  cathode  through  the 
membrane  under  combined  influences  of  water  diffusion  and 
electro-osmotic  drag.  Its  value  is  close  to  zero  as  recently  mea¬ 
sured  by  Liu  et  al.  [44].  In  the  present  work,  for  simplicity  we 
assume  a  equal  to  zero  while  a  non-zero  a  can  be  easily  incor¬ 
porated  in  the  present  model.  Once  the  reaction  rate  is  available, 
the  following  boundary  conditions  can  be  applied  [45]: 

pu-n  =  -(MH2°SU2°  -  M°2S°2)n 

'mff 

Mk 


YcU 


ch2°  . , 


,  „dCf2° 

nLeff  §  I 

g  -  "r 


ycuC°2  ■  n  -  £>°2’eff 


3  n 
dC°2 
3  n 


'"'L  -  |  .  ,1  = 


S°2 


(31) 


4.  Numerical  procedures 


The  governing  equations,  Eqs.  (1),  (2),  (15)  and  (22),  along 
with  their  appropriate  boundary  conditions  are  discretized  by  the 
finite  volume  method  and  solved  by  a  commercial  flow  solver, 
Fluent®  (version  6.0. 12),  using  the  SIMPLE  (semi-implicit  pres¬ 
sure  linked  equation)  algorithm  [46].  The  SIMPLE  algorithm 
updates  the  pressure  and  velocity  fields  from  the  solution  of 
a  pressure  correction  equation,  solved  by  algebraic  multi-grid 
(AMG)  method.  Following  the  solution  of  the  flow  field,  species 
equations  are  solved.  The  source  terms  and  physical  proper¬ 
ties  are  implemented  in  a  UDF  (user-defined  function)  and 
the  species  transport  equations  are  solved  through  the  soft¬ 
ware’s  user-defined  scalars.  The  mesh  of  a  single-channel  PEFC 
employed  here  for  a  numerical  study  is  shown  in  Fig.  8.  Geomet¬ 
rical  and  operating  parameters  are  listed  in  Table  2.  Ten  thousand 
(20  x  100  x  5)  computational  cells  are  used  to  capture  the  com¬ 
plex  two-phase  phenomena  in  the  PEFC  channel.  In  addition, 


Table  2 

Geometrical,  physical  and  operating  parameters 


Quantity 

Value 

Gas  channel  depth 

0.3  mm 

Gas  channel  width 

0.5  mm 

Gas  channel  length 

0.1m 

Cathode  pressure,  P 

1.0  atm 

Average  current  density,  7av 

1.0  A  cm-2 

Temperature  of  fuel  cell,  T 

353  K 

O2  diffusivity  in  the  gas  channel  at  standard 

3.2348  x  10“5  m2  s”1 

condition,  D0;o2 

H2O  diffusivity  in  the  gas  channel  at 

7.35  x  10  5  m2  s  1 

standard  condition,  D0,h2o 

Viscosity  of  liquid  water,  fi \ 

3.5  x  10  4  kg  (ms)  1 

Surface  tension,  liquid-water-air  (80  °C),  0 

0.0625  Nm"1 

Contact  angle  of  the  channel  wall,  6C 

O 

O 

OO 

Liquid-vapor  phase  change  latent  heat,  hfs 

2.26  x  106Jkg-> 

overall  species  balance  is  checked  in  addition  to  the  equation 
residuals  as  important  convergence  criteria.  These  species  bal¬ 
ance  checks  also  ensure  physically  meaningful  results.  In  all  the 
simulations  to  be  presented  in  the  next  section,  values  of  species 
imbalance  are  all  less  than  1%  and  equation  residuals  smaller 
than  10-6. 

5.  Results  and  discussion 

It  is  instructive  to  perform  a  one-dimensional  (along  the  chan¬ 
nel)  analysis  based  on  the  above-described  model.  We  consider  a 
single-straight  channel  and  the  isothermal  condition  with  focus 
on  variations  of  pressure  and  liquid  saturation  along  the  channel. 
The  cathode  side  is  selected  for  analysis  as  shown  in  Fig.  6(b), 
while  the  anode  channel  problem  can  be  derived  similarly. 

5.1.  Buildup  of  liquid  water 

Considering  the  steady  state,  integration  of  Eq.  (1)  in  one- 
dimension  yields: 


AXzpUy  -  (Axz 


f} 

pUy)\in  + 

Jo 


Sm  T^?rnem  d  y 


=  (AxzPUy) |in  +  Lz,me m  ^  7.(!}  dV  MHz 


2  F 


(32) 


610 


Y.  Wang  et  al.  /  Journal  of  Power  Sources  179  (2008)  603-617 


where  ()|in  represents  the  value  of  a  quantity  at  the  inlet  and 
Sm  denotes  the  mass  source  due  to  ORR.  Note  that  the  net  mass 
source  consists  of  water  production  minus  oxygen  consumption. 
Assuming  the  capillary  effect  and  diffusive  transport  along  the 
channel  direction  are  negligible  as  compared  to  the  convection, 
integration  of  Eq.  (15)  for  species  water  results  in: 

Ax-ycityCH2°  =  (AxzyciiyCU2°)\m  +  [  SH2°Lz^memdy 

Jo 

=  (AxzYcUyCU2°)\in  +  (33) 

where  SU2°  denotes  the  water  production  rate  of  ORR.  Substi¬ 
tuting  Eq.  (32)  into  Eq.  (33)  to  cancel  the  mixture  velocity,  uy , 
on  the  left-hand  side  yields: 

ycCH2°  _  (A^MyCH2°)|in  +  Lz,mem((/o'  I(y)dy)/2F) 

P  (Axzpuy)\in  +  I(y)dy)/2F)MU2  (  } 

Using  the  inlet  condition,  Eq.  (27),  one  can  further  express 
the  above  as 


In  the  region  of  Y<  Y0,  there  is  single-phase  flow  in  the  chan¬ 
nel.  Note  that  p  and  pg  represent  the  mixture  density  and  gaseous 
density  respectively  and  their  relation  is  expressed  in  Eq.  (3).  The 
water  concentration  can  be  calculated  through  Eq.  (35)  and  (38). 
Other  variables,  such  as  O2  concentration  and  mixture  velocity, 
can  be  readily  obtained  similarly. 

Assuming  a  constant  current  density,  Y0  can  be  explicitly 
expressed  as 


£c  [(l/C°2)((CsHa|0/pg)p  -  CH2°)]|in 

0  2  1  -  (C^t°/ pg)MH2 


(40) 


As  §c  ->  00,  F0  ->  00  unless  CH2°|in  >  C^f0.  Wang  and 
Wang  [45]  concluded  that  the  electrochemical  reaction  has  little 
effect  (<5%)  on  the  cathode  gas  density  and  flow  velocity  under 
common  operation  of  PEFCs.  Therefore,  Eq.  (39)  can  be  further 
simplified  to  calculate  Y0  by  neglecting  mass  injection  and  gas 
density  variation. 

In  the  region  of  Y>  Y0,  liquid  water  emerges  and  hence  there 
is  two-phase  flow  in  the  channel.  Substituting  Eq.  (38)  into  Eq. 
(35),  one  obtains: 


fe/2)(CH20/C02)|in  +  ((/0F  l(Y)dY)/ /av)  -  (l/Pg)CHa20(fe/2)(p/C02)|in  +  ((/pF  /(F)dF)//av)MH2) 
((1/MH2°)  -  (l/pg)C^°)((£c/2)(p/ C°2)|in  +  ((f0F°  /(f)d?)//av)MH2) 


KcCH2°  =  (fe/2C°2)CH20)|in  +  ((/0y  7(f)d?)//av) 
P  (($c/2C°2)p)\in  +  ((f0y  /(¥)  dP)//av)MH2 


Here,  Y  is  the  dimensionless  distance  from  the  inlet  (see  Fig.  8), 
namely 

Y=  j-  =  (36) 

Ly  ^mem 

The  average  current  density,  /av,  can  be  calculated  by  inte¬ 
grating  the  current  density  profile,  I(Y)\ 

4v  =  I'  I(Y)dY  (37) 

Jo 


Note  that  the  left-hand  side  of  Eq.  (35)  is  solely  determined 
by  liquid  water  saturation  while  the  right-hand  side  is  a  function 
of  operating  parameters,  such  as  §c  and  the  inlet  humidity,  as 
well  as  the  axial  location,  Y.  Given  the  expressions  of  mixture 
density,  Eq.  (3),  and  advection  correction  factor,  Eq.  (16),  the 
left-hand  side  can  be  further  rearranged  as 


KcCH2° 

P 


Ch2° 

Pg 

(U_ 

\M*20 


Y  <Y0 


+  k  +  ~Cf2°,  Y0<Y 


Pg 


(38) 


where  Y0  is  the  transition  point  from  the  single-  to  two-phase 
regions,  i.e.  when  CH2°|y0  =  C^°.  That  is,  Y0  can  be  deter¬ 
mined  by  setting  CH2°  on  the  left-hand  side  of  Eq.  (35)  to  be 
ru  2o. 

^sat  • 

C»2°  =  (fe/2C02)CH2°)|in  +  (( g°  /(F)dF)//av) 

Pg  ((§c/2C°2)p)|in  +  ((/q°  I(Y) d?)//av)MH2 


The  liquid  mobility,  k\,  is  a  function  of  liquid  water  satura¬ 
tion,  s,  as  shown  in  Eq.  (8)  and  Fig.  7.  Through  rearrangement 
from  Eqs.  (8)  and  (9),  one  can  express  liquid  water  saturation 
explicitly  by  the  liquid  mobility,  k\,  as 


S->Sir  _  _ 1 _ 

1  -  %  “  (((1  -  WAzXwgM))17"*  +  1 


(42) 


Substituting  Eq.  (41)  into  Eq.  (42),  we  can  obtain  an  explicit 
solution  of  water  saturation  as  a  function  of  operating  conditions 
and  axial  location.  Note  that  the  effect  of  residual  saturation,  Sjr, 
has  been  accounted  for  in  Eq.  (42). 

Once  liquid  saturation,  s,  becomes  available,  the  mixture  den¬ 
sity  can  be  calculated  through  Eq.  (3),  which  can  further  be 
substituted  into  Eq.  (32)  to  calculate  the  mixture  velocity  by: 

_  (puy)\in  +  ((Jy  I(Y)dY)/2F)Mii2 

Uy  - 

P 

Considering  the  condition  of  a  constant  current  density,  Eq. 
(41)  can  be  simplified  as 


(4/2)(CH2°/ C°2)|in+T— (l/pg)C^°((^c/2)(p/  C°2)|in+TMH2) 
((1/MH20)  _  (l/pg)CsHaf°)(&/2)(p/C°2)|in  +  YM*2) 

(44) 


When  full  humidification  is  applied  at  the  inlet,  Eq.  (44)  can 
be  further  simplified  to: 

x  = _ (1  -  (1  /Pg)C^°MH2) _ 

1  ((1/Mh2)  -  (l/pg)C^°)((fc/2)(p/ C°2)|in  +  FMH 2) 

(45) 

In  the  above  equation,  §c  — >►  00  leads  to  Xc  0  and  hence 
s  — >  as  shown  in  Eq.  (42). 
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For  the  oxygen  equation  of  Eq.  (15),  following  the  similar 
integration,  one  can  arrive  at  the  following  formula: 

AxzycuyC°2  =  Axz(ycuyC02)\in+  T  S°2LZ  ,mem  d y 

Jo 

=  Axz(ycUyC®2)\in  —  fz,mem  .  _  -  (46) 

4  F 

Therefore, 

co2  =  (YcUyC°i)\m-((f*  I(Y)dY)/4F) 

YcUy 

The  right-hand  side  of  the  above  equation  can  be  calculated 
through  the  known  operating  conditions  and  liquid  water  satu¬ 
ration,  .v. 

Liquid  water  saturation,  s,  is  an  important  parameter  char¬ 
acterizing  the  two-phase  flow  in  fuel  cell  channels.  Fig.  9  plots 
the  liquid  water  saturation  profiles  along  the  channel  under  a 
full-humidity  operation  calculated  from  Eqs.  (42)  and  (45)  with 
the  cell  temperature  of  80  °C  and  pressure  of  100  kPa.  We  also 
set  the  irreducible  residual  saturation  s u-  equal  to  zero  to  facili¬ 
tate  explaining  the  results.  The  effects  of  the  irreducible  residual 
saturation  S[Y  can  be  readily  obtained  through  Eq.  (42).  Since  the 
inlet  gas  is  already  fully  humidified,  liquid  water  emerges  at 
the  channel  entrance  due  to  local  water  production.  The  satura¬ 
tion  level  rises  drastically  at  the  initial  stage  followed  by  a  slow 
increase  in  the  downstream.  One  reason  for  this  trend  is  that  the 
mobility  of  liquid  water,  X\,  is  strongly  dependent  on  the  liquid 


Dimensionless  distance  or  Y 

Fig.  9.  Liquid  water  saturation  profiles  along  the  channel,  predicted  from  ana¬ 
lytical  solutions  and  numerical  simulations. 

water  saturation,  as  shown  in  Eqs.  (8)  and  (9),  with  the  value  of 
k\  approximately  proportional  to  snjc  at  low  saturations.  There¬ 
fore,  liquid  mobility  is  very  low  in  low-saturation  regions  and 
liquid  water  builds  up  quickly  as  shown  in  Fig.  9.  In  addition,  it 
can  be  seen  that  the  higher  the  stoichiometric  ratio,  §c,  or  the  gas 
flow  rate,  the  lower  the  saturation  level.  This  can  be  explained 
by  the  fact  that  higher  velocity  of  the  gas  is  efficient  in  draining 
the  liquid  by  overcoming  surface  interactions.  Under  common 
operating  conditions  of  a  fuel  cell,  with  §c  =  2.0,  the  predicted 
saturation  can  reach  as  high  as  20%  towards  the  channel  outlet, 
in  good  agreement  with  experimental  observations.  In  addition, 
3D  simulation  results  of  the  cross-section  averaged  saturation 
are  plotted  in  the  same  figure  and  a  good  agreement  is  obtained 
indicative  of  validity  of  the  assumptions  made  in  derivation  of 
the  analytical  solutions. 


Fig.  10.  Simulation  results  of:  (a)  liquid  water  saturation;  (b)  pressure  distribution  at  the  middle  cross-section  of  the  channel. 


612 


Y.  Wang  et  al.  /  Journal  of  Power  Sources  179  (2008)  603-617 


4e-4  m/s  2e-6  m/s 


Fig.  1 1 .  Simulation  results  of:  (a)  liquid  phase  velocity;  (b)  component  of  liquid 
velocity  in  the  depth  direction  at  the  middle  cross-section  of  the  channel. 


Detailed  simulation  results  are  depicted  in  Figs.  10  and  11, 
displaying  the  distributions  of  liquid  water  saturation,  gas  pres¬ 
sure,  and  liquid  velocity  in  the  mid- section  of  the  channel. 
Fig.  10(a)  indicates  a  small  variation  of  saturation  in  the  depth 
direction.  This  is  due  to  the  small  size  of  the  channel  cross- 
section  compared  with  the  channel  length.  Therefore,  a  slight 
liquid  water  motion  exists  in  the  channel  depth  direction  from  the 
GDL-channel  interface,  where  water  is  injected  into  the  chan¬ 
nel,  to  the  opposing  sidewall  as  shown  in  Fig.  1 1(b).  In  addition, 
Fig.  10(b)  indicates  that  the  gas  pressure  gradient  increases  along 
the  channel  which  can  be  explained  by  increase  of  liquid  satura¬ 
tion  along  the  channel  as  shown  in  Fig.  10(a).  Fig.  11(a)  shows 
that  the  liquid  flow  along  the  channel  speeds  up,  due  to  con¬ 
tinual  water  injection  from  the  GDL-channel  interface  (the  left 
wall)  due  to  ORR.  As  the  gas  phase  is  already  fully  saturated  at 
the  inlet,  the  fuel  cell  channel  totally  relies  on  the  liquid  phase 
to  remove  product  water.  In  addition,  liquid  axial  velocity  is 
small  along  the  channel  with  the  magnitude  of  10-4  ms-1.  This 
velocity  magnitude  renders  a  time  constant  ~  1000  s  for  water 
drainage  or  the  reverse  process  of  liquid  water  accumulation  in 
the  channel,  given  the  channel  length  of  0. 1  m.  This  residence 
time  is  consistent  with  experimental  observations  on  fuel  cell 
transients  [47]  and  much  larger  than  the  ones  for  gas  diffusion, 
membrane  hydration  [48-50],  GDL  drying  [17,5 1]  and  cold  start 
[52,53]. 


Fig.  12.  Analytical  results  of  liquid  water  saturation  profiles  under  non-uniform 
current  density  distribution. 


Fig.  12  shows  liquid  saturation  profiles  under  varying  current 
density  distributions  as  calculated  from  the  analytical  solution, 
Eqs.  (41)  and  (42).  As  current  density  is  proportional  to  water 
production  rate,  the  current  density  distribution  has  profound 
impacts  on  the  local  phase  distribution  as  shown  in  Fig.  12. 
However,  the  average  current  densities  are  set  the  same,  therefore 
the  same  level  of  liquid  water  saturation  is  predicted  at  the  outlet 
for  all  the  cases  despite  of  different  local  current  density. 

While  the  above  results  focus  on  fully  saturated  gases  at  the 
inlet,  low-humidity  operation  is  more  common  in  the  automotive 
application  of  PEFC.  Fig.  13  shows  simulated  profiles  of  liquid 
water  saturation  under  different  inlet  stoichiometric  ratios  for 
the  inlet  gas  dew  point  of  70  °C.  Because  low-humidity  gas  is 
fed  in,  there  is  no  liquid  water  in  the  channel  entrance  region 
and  channel  flow  is  single-phase  in  nature.  As  the  channel  con¬ 
tinues  receiving  water  produced  from  ORR  via  the  GDL,  the 
water  concentration  in  the  gas  increases  till  saturation.  After 
that,  water  addition  leads  to  condensation  as  the  gas  is  already 
saturated.  This  is  the  transition  point  from  single-  to  two-phase 
flow  in  channels,  as  clearly  shown  in  Fig.  13.  At  the  initial 
stage  of  two-phase  flow,  the  liquid  saturation  increases  gradually 
due  to  the  fact  that  a  portion  of  channel  flow,  near  the  sidewall 
opposing  to  the  GDL,  is  in  pure  single-phase  region  which  holds 
water  to  increase  its  water  content.  After  the  entire  channel  is  in 
the  two-phase  regime,  the  liquid  saturation  experiences  a  sharp 
rise  followed  by  a  slow  increase  towards  the  channel  outlet.  In 
addition,  as  the  air  stoichiometry  increases,  the  transition  from 
single-  to  two-phase  flow  delays,  due  to  the  larger  capacity  of 
gas  to  uptake  water  vapor  before  reaching  saturation.  This  con- 


Fig.  13.  Simulation  results  of  liquid  water  saturation  distributions  under  low- 
humidity  inlet  conditions. 


Y.  Wang  et  al.  /  Journal  of  Power  Sources  179  (2008)  603-617 


613 


Fig.  14.  Simulation  results  of:  (a)  liquid  velocity;  (b)  component  of  liquid  veloc¬ 
ity  in  the  depth  direction;  (c)  component  of  gas  velocity  in  the  depth  direction 
at  the  middle  of  the  gas  channel. 


elusion  can  also  be  drawn  from  the  analytical  solution  for  Y0, 
Eq.  (39). 

Fig.  14  shows  the  simulation  results  of  liquid  velocity,  veloc¬ 
ity  components  of  liquid  and  gas  flow  in  the  depth  or  x-direction 
for  the  inlet  stoichiometry  of  2  and  dew  point  temperature  of 
70  °C.  It  can  be  seen  from  Fig.  14(a)  that  there  exists  a  region 
where  single-  and  two-phase  flows  coexist  and  the  liquid  flow 
accelerates  after  the  entire  channel  is  in  the  two-phase  regime. 
Fig.  14(b)  clearly  demonstrates  that  the  present  M 2  model  can 
resolve  differing  velocity  fields  of  gas  and  liquid,  indicating  that 
it  is  a  full  two-fluid  model.  The  liquid  flow  in  the  channel  depth 
direction  is  away  from  the  GDF  surface  while  the  gas  flow  is 
towards  the  GDF,  as  shown  in  Fig.  14(c).  The  magnitude  of  the 
transverse  gas  velocity  is  small  (~10-3ms-1)  compared  to  the 
axial  gas  velocity. 

Fig.  15  displays  the  liquid  saturation  profiles  for  different  dew 
point  temperatures  of  inlet  gases  at  the  stoichiometry  of  2.  It  can 
be  seen  that  the  higher  the  gas  dew  point,  the  earlier  transition 
into  the  two-phase  flow,  as  expected. 


0 


0 


Fig.  16.  Schematic  of  branch-merge  flow  field  and  the  corresponding  porous 
channel  model  (K\  >Kf).  The  cycle  shadow  represents  the  liquid  water  trapping 
at  the  heterogeneity. 


5.2.  Two-phase  flow  in  heterogeneous  channels 


The  foregoing  results  and  discussion  have  been  focused  on 
channels  of  uniform  cross-section  and  with  constant  surface 
wettability,  while  many  PEFC  designs  have  non-uniform  chan¬ 
nels  which  feature  either  geometrical  heterogeneity  or  surface 
wettability  variation  along  the  flow,  or  both.  The  two  types  of 
heterogeneity  are  considered  together  and  regarded  generally  as 
a  spatial  variation  in  channel  permeability  (or  hydraulic  conduc¬ 
tance)  and  surface  contact  angle.  For  example,  a  branch-merge 
channel  can  be  treated  as  a  porous  medium  with  differing  per¬ 
meability  or  hydraulic  conductance,  as  shown  in  Fig.  16.  This 
branch-merge  flowfield  is  a  popular  design  in  fuel  cells  due  to  its 
ability  to  mix  and  re-distribute  reactants  in  channels  thereby  alle¬ 
viating  a  tendency  of  flow  maldistribution.  However,  the  channel 
design  may  also  lead  to  water  trapping  at  the  transition  point  as 
shown  in  Fig.  16.  In  addition,  a  branch  of  parallel  channels  in  a 
fuel  cell  usually  share  the  same  inlet  or  outlet,  also  called  mani¬ 
fold  as  shown  in  Figs.  3  and  4,  which  resembles  either  the  merge 
or  branched  pattern. 

For  liquid  water  transport  through  a  homogeneous  porous 
channel,  the  capillary  action  can  be  written,  through  substituting 
Eq.  (19)  into  the  second  term  on  the  right-hand  side  of  Eq.  (15), 
as 


Al  Act 

~^KVPC 

V 

(48) 


Fig.  15.  Analytical  solutions  of  liquid  water  saturation  profiles  under  low- 
humidity  inlet  conditions. 


The  above  term  can  be  treated  as  a  diffusion  term  with  the 
following  form: 


V  • 


\  Mk  pg) 


j\ 


pg  /  v  dv 


(49) 


Thus,  the  capillary  force  drives  water  from  a  higher  saturation 
to  a  lower  saturation  region  in  a  homogenous  porous  medium. 
However,  in  a  heterogeneous  porous  medium,  Eq.  (48)  must  be 
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02+4e+4H+^2H20+heat 


<=> 

Air+H20 


l=> 

Air+H20 


Fig.  17.  Schematic  of  permeability  heterogeneity  in  fuel  cell  channels. 


Dimensionless  distance  or  Y 

Fig.  18.  Simulation  result  of  the  liquid  water  saturation  around  the  heterogeneity 
shown  in  Fig.  17. 


rewritten  in  the  following  general  form: 


V  • 


Ca\  A-lA-g 

——^KVPc(cr,  6C,  s,  K,  s) 
p%)  v 


=  V  • 


AqA.p 


dPc 

da 


K  (  ^V0C 


dPo 

dOc 


dPc  dPc  dPc 

+  —  Vs  +  —  VK  +  —Vs 
ds  dK  ds 


(50) 


It  can  be  seen  that  spatial  variations  in  surface  tension  a,  con¬ 
tact  angle  0C,  porosity  s  and  permeability  K  add  extra  terms  into 
the  capillary  action  in  addition  to  the  capillary  diffusion  under 
saturation  gradient.  For  simplicity,  here  we  only  consider  the 
heterogeneity  in  the  permeability  (or  hydraulic  conductance), 
K ,  while  keeping  other  properties  constant.  Then  the  additional 
term  due  to  geometrical  heterogeneity  can  be  expressed  as 


V  • 


v 


=  V  • 


XlV 1/2 

V 

(51) 


Further  assuming  a  linear  profile  of  K  along  the  channel,  as 
shown  in  Fig.  17 


K  =  k\y  +  k2 


(52) 


This  localized  phenomenon  can  be  readily  explained  by  Eq.  (50) 
where  the  extra  term  acts  like  an  additional  capillary  force  to  hold 
up  water.  On  the  global  scale,  since  the  extra  terms  in  Eq.  (50) 
are  mathematically  in  conservative  form,  there  is  no  extra  water 
added  besides  water  production  by  ORR.  Finally,  it  is  notewor¬ 
thy  that  if  reaching  a  high-liquid  saturation,  the  localized  spike 
may  cause  channel  clogging  by  liquid  water  and  flow  shutdown 
in  a  multiple,  parallel  channel  configuration. 

It  should  be  noted  that  while  Eq.  (50)  is  derived  for 
isotropic  porous  media,  a  similar  analysis  could  be  performed  for 
anisotropic  media.  More  discussion  on  two-phase  flows  in  het¬ 
erogeneous  porous  media  can  be  found  in  the  field  of  petroleum 
and  geological  systems  [26,32]. 

5.3.  Two-phase  pressure  drop 

As  the  presence  of  liquid  water  hampers  gas  flow  in  fuel  cell 
channels,  the  gas  pressure  gradient  in  two-phase  flow  is  greater 
than  that  in  single-phase.  The  increased  pressure  drop  due  to  liq¬ 
uid  water  accumulation  is  a  key  cause  for  flow  maldistribution, 
as  schematically  shown  in  Fig.  19.  Liquid  water  in  the  lower 
channel  incurs  additional  pressure  drop  which  forces  the  gas 
flowrate  to  be  reduced  in  the  lower  channel  and  increased  in 
the  upper  channel  in  order  to  maintain  the  same  pressure  drop 
between  the  inlet  and  outlet  manifolds.  Such  a  two-phase  flow 
maldistribution  dramatically  reduces  PEM  fuel  cell  performance 
as  well  as  durability. 


One  can  reach 


KV(K~l/2)  =  K 


f0  \ 

-b^-(3/2) 

2 

Vo  / 


(°  \ 

£-0/2) 

\0  / 


(53) 


(a) 


Fig.  18  shows  the  liquid  water  saturation  profile  along  a  het¬ 
erogeneous  channel  schematically  shown  in  Fig.  17.  It  can  be 
seen  that  the  saturation  profile  remains  almost  unchanged  in  the 
uniform  sections  of  the  channel  while  a  spike  is  predicted  in  the 
heterogeneous  region.  This  result  indicates  that  the  influence 
of  local  heterogeneity  affects  only  the  local  phase  distribution. 


tr  <b> 

Fig.  19.  Schematic  of  flow  in  multiple  channels:  (a)  uniform  distribution;  (b) 
maldistribution.  The  black  shadowed  regions  in  (b)  denote  the  liquid  phase. 
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Fig.  20.  Two-phase  pressure  factor  vs.  the  stoichiometric  flow  ratio  under  dif¬ 
ferent  inlet  humilities  for  a  single  channel  ( n^  =  4). 
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Fig.  21.  Experimental  validation  of  the  channel  two-phase  model. 


The  gas  pressure  gradient  is  linked  to  the  gas  velocity  via  the 
Darcy’s  law,  namely 


K(  1  -  s)nk  dpg 
Mg Ly  d Y 


(54) 


In  view  of  the  nitrogen  content  of  ~79%  in  air,  we  can  safely 
assume  that  the  gas  velocity  remains  nearly  constant  along  the 
channel,  equal  to  the  inlet  value,  Mg  ~  .  Indeed,  Wang  and 

Wang  [45]  numerically  showed  that  variation  of  the  gas  phase 
velocity  is  small  (<5%)  in  common  PEFC  operation.  Thus,  Eq. 
(54)  can  be  rearranged  to: 


dpg  _  fav^mem  1  Mg^y  £c 

~dY  ~  AF  C°2 Axz  K  (l-s)nk 
If  a  dimensionless  pressure  is  defined  as 


(55) 


|AP|  = 


_ |Agg| _ 

(/ayAmem  /AF)(\ / C°2  Axz)(n„Ly/ K) 


(56) 


One  obtains  that 

I  AP|  =  §c 


r _ _ 

Jo  (1  -  S(Y,Zc 


RHi 


d  Y  =  |  A />(£c)| 


(57) 


Eq.  (57)  shows  that  the  dimensionless  pressure  drop  is  only 
a  function  of  the  air  stoichiometry  §c  and  inlet  humidity,  RH.  In 
the  single-phase  regime,  Eq.  (57)  is  reduced  to: 

|  AP1^!  =  £c  (58) 


Thus,  a  two-phase  pressure  factor,  <f>,  can  be  defined  as 


|AP| 

\AP1<t>\ 


f1  1 

/  - dr 

Jo  [1  -  s(Y,  $C)  RHin)]"‘ 


(59) 


Eq.  (59)  indicates  that  the  two-phase  pressure  factor  depends 
only  on  the  air  stoichiometry  and  inlet  RH.  Fig.  20  shows  the 
analytical  result  of  the  two-phase  pressure  factor  as  a  function 
of  stoichiometric  ratio  for  different  inlet  dew  points  or  RHs.  It 
can  be  seen  that  at  low- stoichiometric  ratios  or  low-channel  flow 
rates,  the  two-phase  pressure  factor  is  greater  than  unity  due  to 


liquid  water  presence.  As  the  flow  rate  increases,  the  occurrence 
of  two-phase  flow  shifts  more  towards  the  downstream  and  hence 
the  two-phase  pressure  factor  is  lower  and  approaches  unity. 
Also,  Eq.  (59)  clearly  indicates  that  the  exponent,  ng,  is  a  key 
factor  determining  the  two-phase  pressure  factor.  Therefore,  one 
can  measure  the  two-phase  pressure  drop  across  the  fuel  cell 
channels  in  order  to  calibrate  ng,  as  to  be  explained  briefly  in  the 
next  sub- section. 

5.4.  Experimental  validation 

To  experimentally  validate  the  above-described  channel  two- 
phase  model  and  to  better  understand  the  applicable  range  of  ng 
specifically  for  fuel  cell  channels,  a  seven-parallel-channel  fuel 
cell  of  0. 1  m  long  was  designed  and  pressure  drop  measurements 
were  carried  out.  While  experimental  details  have  been  docu¬ 
mented  by  Hussaini  and  Wang  [55]  and  thus  are  not  repeated 
here,  the  experimental  data  of  two-phase  pressure  drop,  selected 
over  a  range  of  air  stoichiometry  and  current  density  to  be  rep¬ 
resentative,  are  plotted  in  Fig.  21,  after  excluding  any  influence 
from  inlet  and  outlet  manifolds.  These  average  data  for  straight 
channels  can  then  be  compared  with  the  present  model  prediction 
for  a  single  channel,  as  shown  in  Fig.  21.  It  is  seen  that  there  is 
slight  under-prediction  using  ng  =  4.  However,  the  experimental- 
model  agreement  is  good  for  both  ng  =  4.5  and  5.0.  It  can  thus  be 
concluded  that  the  present  channel  two-phase  model  is  reason¬ 
ably  accurate  and  the  exponent  in  the  liquid  permeability  likely 
lies  between  4.5  and  5.0. 

6.  Conclusions 

In  this  work,  we  have  developed  a  two-phase  model  for 
channel  flow  in  a  PEM  fuel  cell,  based  on  the  M 2  formalism.  The¬ 
oretical  analyses  were  performed  to  calculate  the  liquid  water 
saturation,  onset  location  of  two-phase  flow,  and  species  concen¬ 
trations  along  the  channel.  Dependence  of  these  physical  profiles 
vital  to  cell  performance  on  key  operating  parameters  such  as  air 
stoichiometry  and  relative  humidity  is  explored.  An  analytical 
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solution  was  developed  to  express  the  liquid  water  saturation 
as  a  function  of  the  stoichiometric  ratio,  relative  humility,  and 
axial  location.  The  model  was  further  implemented  for  three- 
dimensional  numerical  simulations  of  two-phase  flow  in  a  single 
channel.  The  results  reveal  that  liquid  water  builds  up  quickly  at 
the  entrance  region  followed  by  a  slow  increase  downstream, 
under  full-humidification  inlet  conditions.  The  predicted  liq¬ 
uid  saturation  reaches  as  high  as  20%.  Subsequently,  two-phase 
flow  in  heterogeneous  channels  were  theoretically  analyzed  and 
numerically  simulated.  Water  trapping  around  the  geometri¬ 
cal  heterogeneity  was  found.  In  addition,  a  two-phase  pressure 
factor  was  defined  which  could  help  establish  a  fundamental 
understanding  of  two-phase  flow  maldistribution  in  fuel  cell 
channels.  Finally,  preliminary  validation  of  the  present  channel 
two-phase  model  against  experimental  pressure  drop  data  was 
carried  out,  with  good  agreement.  Future  work  is  to  validate  the 
present  model  against  experimental  liquid  water  saturation  pro¬ 
files  once  they  are  available,  as  well  as  to  integrate  the  channel 
two-phase  model  into  a  full  cell  model. 

Appendix  A.  Nomenclature 


A  area  (m2) 

Ck  molar  concentration  of  species  k  (mol  m-3) 

D  mass  diffusivity  of  species  (m2  s-1) 

F  Faraday’s  constant  (96,487  C/equivalent) 

I  current  density  (A  cm-2) 

k  thermal  conductivity  (W  (m  K)  - 1 ) 

kv  relative  permeability 

K  permeability  (m2) 

L  length  (m) 

mfk  mass  fraction  of  species  k  in  liquid  phase 

M  molecular  weight  (kg  mol- 1 ) 

n  the  direction  normal  to  the  surface 

P  pressure  (Pa) 

R  gas  constant  (8.134  J  (mol  K)-1) 
s  stoichiometry  coefficient  in  electrochemical  reaction 

or  liquid  saturation 

S  source  term  in  transport  equations 

t  time  (s) 

T  temperature  (K) 

u  velocity  vector  (ms-1) 

X  mole  fraction 

Greek  letters 

a  net  water  transport  coefficient  per  proton 

£  porosity 

k  ionic  conductivity  (S  m- 1 ) 

Xk  mobility  of  phase  k 

v  kinematic  viscosity  (m2  s-1) 

§  stoichiometric  flow  ratio 

p  density  (kg  m- 3 ) 

a  surface  tension  (N  m- 1 ) 

r  tortuosity  factor 

0  two-phase  pressure  factor 


Superscripts  and  subscripts 
a  anode 

av  average 

c  cathode 

eff  effective  value 

g  gas  phase 

in  inlet 

k  species 

1  liquid  phase 

mem  membrane 

o  standard  condition,  273.15  K  and  101.3  kPa  (1  atm) 

sat  saturate  value 
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